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ABSTRACT 

We study the growth of black holes (BHs) in galaxies using three-dimensional smoothed particle 
hydrodynamic simulations with new implementations of the momentum mechanical feedback, and 
restriction of accreted elements to those that are gravitationally bound to the BH. We also include 
the feedback from the X-ray radiation emitted by the BH, which heats the surrounding gas in the 
host galaxies, and adds radial momentum to the fluid. We perform simulations of isolated galaxies 
and merging galaxies and test various feedback models with the new treatment of the Bondi radius 
criterion. We find that overall the BH growth is similar to what has been obtained by earlier workers 
using the Springel, Di Matteo, & Hernquist algorithms. However, the outflowing wind velocities and 
mechanical energy emitted by winds are considerably higher (v w ~ 1000 — 3000 km s~ x ) compared 
to the standard thermal feedback model (v w ~ 50 — 100 km s~ x ). While the thermal feedback model 
emits only 0.1% of BH released energy in winds, the momentum feedback model emits more than 30% 
of the total energy released by the BH in winds. In the momentum feedback model, the degree of 
fluctuation in both radiant and wind output is considerably larger than in the standard treatments. 
We check that the new model of the BH mass accretion agrees with analytic results for the standard 
Bondi problem. 

Subject headings: accretion, accretion disks - black hole physics - galaxies: active- galaxies: nuclei - 
galaxies: starburst - quasars: general 



1. INTRODUCTION 

Accretion of matter onto a supermassive black hole 
(SMBH) is believed to emit enormous amounts of elec- 
tromagnetic, luminous radiation a nd drive powerful 
jets, winds or outflows of quasar (iLvnden- Bcll 1969; 
Rees 1984). The energy outputs emitted by accret- 
ing SM BHs at the centers of bulges and elliptical 
galax ies (jKormendv fc Richstond 119951 : iRichstone et al.1 
[1998Q are believed to play an important role during 
galaxy formation and evolution, as revealed by many 
empirical correlations between their masses Mbh and 
host galaxy properties, e.g . , the dispersion a o f the 
host galaxy (iGebhardt et alj|2000t iTremaine et alj|2002t 
IGultekin et all 120091: iGraham et all 1201 lh. bulge stellar 
mass dDresslerl 119891: iKormendvl 119931: iMagorrian et al.l 
[1991 iMarconi fc H untl 120031: I Gultekin et al l 120091) and 
the bulge binding energ y (|Aller fc RichstoneT l2007t 
Barway & Kcmbhavi 2007). These correlations have led 
to the development of models where the SMBHs are 
linked by feedback from the central SMBHs, i.e., feed- 
back via the mass ejection by winds or jets, or the emit- 
ted radiation regulates the mass accretion rate and the 
final SMBH mass. The feedback can be i n the forms of 
radiatively or me c hanic ally driven winds (|FabianlH99l 
iProga et al.l 120001 [2008), or of radiative effects such as 
Compton heating and photo ionization (Sa zonov et al.l 
120051: iCiotti fc Ost riker 2007), or of radiation pressure 
(jDeBuhr et all 120101 l2011aT i. Thermal feedback, heat- 
ing by an unspecified mechan ism has also been em - 
ployed by many authors, e.g., Springcl et al.l (|2005b|) : 



iHopkins et all (|2005f ): I Johansson et al.l (j2009al |bh. 

In this rapidly developing field of active galactic nu- 
cleus (AGN) feedback, it is natural for large scale nu- 
merical simulations to explore first the simplest schemes. 
Thus physical processes, which cannot all be included 
in the first three-dimensional hydrodynamic simulations, 
have b een treated in a select ive fashion. For exam- 
ple, in ISpringel et all (I2005bl here after "SdMH05"), 
among the first who reproduced the AGN feedback in 
a three-dimensional smoothed particle hydrodynamics 
(SPH) code, the feedback was assumed to be purely ther- 
mal. That is, some fraction of the bolometric luminos- 
ity of accreting BHs is deposited as thermal energy to 
the neighboring gas particles via a mechanism that is 
not specified. The authors found that this thermal feed- 
back treatment regulates and then terminates the fur- 
ther growth of the BH and expels gas from the central 
region in the galaxy in a quasar-driven outburst. This 
pioneering work, however, did not specify how the in- 
jected energy reaches the thermally heated gas particles. 
The conveyance of the energy is likely to be via cither 
radiation or a wind, and in both cases momentum must 
be added with the thermal energy. If the energy is trans- 
ferred via a wind with velocity i> w , then the mass transfer 
may be significant and, since the ratio of momentum to 
energy scales as l/u w , the momentum transfer is corre- 
spondingly larger. 

In a first attempt to study the relative impor- 
tance of the different processes, utilizing one- and 
two-dimensional computations, it was found that mo- 
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mentum injection is the dominant mode of feedback 
HOstriker et al.l l2010l here after "OCCNP10"), because 
the very short cooling time of dense gas makes thermal 
input relatively inefficient. At the high densities of the 
central region of the galaxies, the cooling time of gas is 
sufficiently short in high resolution simulations, so that it 
cannot retain the injected thermal energy and efficiently 
convert it to kinetic energy. However, if the input to 
the surrounding fluid is via winds, then the return of 
mass and momentum to that fluid can be the dominant 
drivers which can reduce the accreted SMBH mass by 
up to a factor of 100 (OCCNP10). Other recent work by 
iDeBuhr et all (|2010l l2011al) . emphasizing the momen- 
tum input from optically thick radiation fields (~ tL/c 
with r ~ 10) has also found the dramatic effects of mo- 
mentum input. 

Observationally, aro u nd 15%-20% of b right quasars 
(iHewitt fc Foltd l200l iDai et al.l 120081 iGibson et all 
2009) show outflows, where blueshifted broad absorp- 
tion lines (BAL) are attributed to subrelativistic (~ 
10,000 km s _1 ) mass ejection. According to the frac- 
tion of quasars with BAL winds together with the as- 
sumption that all quasars have outflows, the BAL out- 
flows are wide, covering at least 20% of the solid an- 
gle. The radia tion-driven winds f r om th e detailed hydro 
simulation by iProga fc KallmarJ (|2004D cover ~ n str, 
supporting this observational estimate. These outflows 
inject mass, momentum, and energy into the surround- 
ing gas and are believed to be more efficient feedback 
agents of the host galaxies than relativistic jets, which 
drill through the ambient gas and put most of the en- 
ergy into the intergalactic medium. 

However, the effects of outflows to the host galax- 
ies and larger cosmological structures are just beginning 
to be studied, because their mass ejection rates, ener- 
getics and sizes were largely undetermined. Adopting 
the recent measurements of mass ejection rates, and ki- 
netic luminosities of the outflows from the absorption- 
line observations (|Arav et al.1 120081: iMoe et al.l 120091 : 
iDunn et alj |2010), many authors introduced the kinetic 
AGN fe edback in adaptive mesh refinement (AMR) simu- 
lation s (|Omma et all2004lKim et al.l201lUDubois et all 
20091) and in SPH si mulations (jNavakshin fc Powerl l2010l: 



DcBuhrelalj[2012). 



It has been shown that the average spectrum of 
AGN has a strong secondary peak in the high-energy 
X-rays (~ 50-100 keV) which is the main contribu- 
tor t o the Compton and secondary metal line heat- 
ing (jSazonov et al.1 120041 ). Some recent studies us- 
ing one- and two-dim e nsional simulation s including 
iCiotti fe Ostrikerl (|2007D : iNovak et all (|2011| ) considered 
the effect of momentum input, heating and radiation 
pr essure from the AGN r adiation. The recent SPH work 
bv lHambrick et aTl ()201lD found that X-ray feedback, the 
heating and radiation pressure associated with the X-ray 
radiation field emitted from the BH, is more effective at 
suppressing star formation and BH mass growth com- 
pared to the traditional thermal feedback approach. The 
X-ray he ating was also fou nd to be effective in AMR sim- 
ulations (|Kim et al.l 1201 lh ■ as it heats the surrounding 
interstellar medium (ISM), keeps it hot for an extended 
duration of time, and effectively self-regulates the growth 
of the BH. 



In many previous numerical studies of galaxy forma- 
tion and evolution that only can resolve hundreds of pc 
to kpc scales, the accretion of gas onto BHs at the cen- 
ters of AGNs occurs on unres olved scales. The Bondi- 
Hoyle-Lyttleton formalism (|Hovle fc Lyttletonl 119391 : 
iBondi fc HovldfTMl : lBondilll952h is commonly adopted, 
to obtain the BH mass accretion rate from the resolved 
larger scale properties of the ambient gas. This assump- 
tion has been made in studies of isolated galaxy and 
merging galaxies (e.g SdMH05, lYounger et all 120081 : 
Uohansson et al. | l2009bf) and as well in cosmological sim- 
ulations (e.g..lSiiacki et all 120071: iDi Matteo et all 1200 
Booth fe S chavc 2 0091 : iTevssier et alll2011l : lDubois et al 



20121) . 



However, in most of the cases treated to date, the SPH 
smoothing length or the numerical resolution of the re- 
spective code is larger than the Bondi radius, i.e., the 
inner flow is numerically unresolved. The limitation of 
numerical and spatial resolution may introduce the un- 
usual and unwanted effect — the neighboring gas parti- 
cles around a BH at any distance from the BH can be 
accreted even though they are not gravitationally bound 
to the BH. This physically awkward treatment of the BH 
mass accretion may result in the very large growth of the 
central BH if we used the standard accretion algorithm 
without feedback. We will introduce a new treatment 
of gas particle accretion, "Bondi radius criterion" , which 
statistically limits the accretion of mass to the gas par- 
ticles which are within the Bondi radius. 

The purpose of the current paper is to introduce (1) a 
modeling of AGN mechanical feedback via winds as ob- 
served in BAL systems that includes mass and momen- 
tum feedback, as well as thermal input; (2) the detailed 
treatment of radiative effects, i.e., X-ray radiative feed- 
back; and (3) a modified BH accretion rate prescription 
using the Bondi radius criterion in the parallel TreeSPH- 
code GADGET. In this paper, we restrict our exploration 
to the simulations of isolated disk galaxies and merging 
galaxies, in order to better understand the specific prop- 
erties of the new treatments. We reserve the follow-up 
papers for detailed analysis of merger simulations and 
cosmological simulations, to be compared with the re- 
sults of observational papers. 

This paper is structured as follows. In Section [21 we 
describe the simulation code, and discuss the BH feed- 
back and new BH accretion prescription. We present 
the results and comparisons between our new model and 
standard prescription in Section [3l We summarize and 
discuss our findings in Section 2] 

2. THE MODELS 

2.1. Numerical Code 

We perform the simu lations using the parallel 
TreeSPH-code GADGET (lSpringelll2005l). T he code em- 
ploys the Lagrangian SPH (see iMonaghanl [T992) tech- 
nique and solves the equations of motion for the col- 
lisionless dark matter and star particles and gas. We 
include the radiative coo ling for a primor dial mixture 
of hydrogen and helium (jKatz et all 119961 ) and a spa- 
tially uniform time-inde pendent local UV background 
(jHaardt fe Madaul 119961) . The gas of the ISM is as- 
sum ed to be a two-phase medium of hot and cold 
gas (jMckee fe Ostrikerlll977HSpringel fe Hernquistll2003l) 



Radiative and Mechanical AGN Feedback 



3 



and stars form from a cold component embedded in suf- 
ficiently dense gas, i.e., p > p t h with the short-lived stars 
supplying an energy of 10 51 erg to the surrounding gas 
per supernovae (SNe). This SN feedback heats the hot 
phase of the ISM and evaporates cold clouds, establish- 
ing a self-regulation cycle for star formation. SN-driven 
galactic winds are not included in this study. 

We include all the basic aspects of the model for black 
hole (BH) acc r etion and feedback adopted in SdMHOS, 
iSpringel etafl (|2005a[HJohansson et al.1 (|2009b[ ). and im- 
plement the momentum and mass feedback. 

2.2. Black Hole Feedback Model 

In the widely adopted SdMH05 model, the numeri- 
cally unresolved accretion onto the BH is related to the 
large scale resolved gas distrib ution using a Bondi-Hoyle- 
Lyttleton parameter i zation (IHovle fe Lvttletonl 119391 : 
IBondi fc Hovlel H941 IBondil [19521 7" In very high res- 
olution treatments where the Bondi radius (Rb ~ 
2GM B h/(c 2 + v 2 ) 1 / 2 ) is r e solved such as ICiotti et al.l 
(pM^ : lNovaket all pOll : iBarai et al.l (|2011j), there is 
no need to specify the accretion rate as the hydro code 
with an appropriate inner boundary condition will cor- 
rectly calculate the accretion rates. The accretion rate 
onto the BH in unresolved flows is estimated as 



Mr = 



47raG 2 M| H p 
(c 2 + v 2 ) 3 / 2 ' 



(1) 



where p and c s are the density and sound speed of the 
surrounding gas, respectively, v is the velocity of the BH 
relative to the surrounding gas. Here a is a dimension- 
less parameter setting the efficiency of the accretion and 
it is conventionally set to be a = 100 in SPH simula- 
tions on the grounds that the low spatial resolution cur- 
rently available would otherwise limit the accretion rate 
to lower than the true value. Adopting a = 100 gives 
reasonable results for t he low resolution s imulati ons dis- 
cussed in SdMH05, and Uohansson et "all (j2009b[ ). but in 
general a is resolution dependent. We adopt different 
values of a for the different resolutions, and details will 
be presented later in Section [3TTI Note that Equation (JTJ) 
describes the accretion onto a point mass surrounded by 
adiabatic (7 = 5/3) gas with properties p, c, and v (in 
Equation (JXJ) ) far away from the BH (r —> 00). Usually 
it is also assumed that the maximum accretion is limited 
to the Eddington rate given by 



M 



(id 



47tGMbhto p 



e r <7 T c 



(2) 



Here m p is the proton mass, <tt is the Thomson cross- 
section and e r is the radiative efficiency assumed to be a 
fixed value of 0.1 adopted from the mean value for radia- 
tively efficient fShakura fc Sunyaevl (pJ373) accretion onto 
a Schwarzschild BH. This value is also sup ported by the 
glob al mass-energy r e lation pointed out bv lSoltanl (|1982j ) 
and lYu fe Tremaind (|2002ft . The accretion rate is then 

A/i n f = min(MB, M cc id)- The BHs are represented by col- 
lisionless "sink" particles, which feel only gravitational 
forces. The properties, including density, temperature, 
and the bulk velocity of the local gas around the BH, 
which then define the accretion rate, are estimated in a 
similar fashion to normal SPH particles. 



Technically, the actual accretion of gas particles onto 
the BH particle is implemented using a stochastic ap- 
proach (SdMH05). For each gas particle j around a BH, 
the probability of being absorbed by the BH is calculated 

as 

WjM ini At 

Pj = - J —z . ( 3 ) 

where Wj is the kernel weight of the gas particle relative 

to the BH, Mint is the BH accretion rate and, At is the 
time step. The gas density p is measured at the position 
of the BH as 

JV 



rriiWi, 



(4) 



where mi denotes the gas particle mass, and N is the 
number of neighboring particles to the BH. From Equa- 
tions ((31) and ([J}, we see that the probability of the jth 
particle being accreted in time interval At is close to 
Pj = Mi n f At/wgas, but influenced by the smoothing ker- 
nel of gas particle near the BH. The gas particle is swal- 
lowed by the BH when the probability pj is larger than 
the generated random number uniformly distributed in 
the interval [0,1]. 

2.2.1. Thermal Feedback 

In previous GADGET based studies, the feedback en- 
ergy from the BH _Ef cc d has typically been assumed to be 
some fraction et of the radiated luminosity L r and cou- 
ples thermally and isotropically to the surrounding gas 

as, 



-Efecd — CfL T = e/e r Mi n fC . 



(5) 



Many authors (e.g., SdMH05) adopt a fixed value of 
ef = 0.05 so that e w = e/e r = 5 x 10~ 3 , i.e., 0.5 per- 
cent of the accreted rest mass energy in total is avail- 
able as thermal energy feedback. Utilizing this value for 
thermal feedback and adopting a — 10 0, the normaliza- 
tion o f the Mbh — cr was recovered by iDi Matteo et al.l 
p005h in disk mergers simulations having a spatial res- 
olution of e gas = 0.1ft. -1 kpc and a mass resolution of 
M ga8 ~6x 10 5 Mq. 

The corresponding feedback energy is distributed as 
thermal energy to the surrounding ~ 64 gas particles 
weighted by the SPH kernel. Note that the results may 
depend on the details of the numerical implementation. 
If ~ 128 particles are used of a given mass rather than 
~ 64, the sound speed c 2 would have been lower with 
corresponding increase of the accretion rate Mb in Equa- 
tion ([TJ. What matters is the mass of the material into 
which the thermal energy is dumped. Thus, if the mass 
per particle were halved and the particle number was 
doubled the results would not change. In this standard 
approach, neither mass nor momentum is added to the 
ambient fluid by the BH and all energy that is added is 
via thermal energy. 

2.2.2. Momentum Feedback 

Accreting BHs are observed to emit broad-line winds 
that convey mas s, energy, and momentum into the 
surrounding gas (tdeJCoolTeUaL| 120011: iMoe et al.l 120091: 
Dun n et all [2010h and our goal is to include these ob- 
served flows in our numerical treatment. To implement 
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mechanical momentum and mass feedback, we define the 
inflowing and outflowing mass rates to be (M- m {, M out f), 
and we use the following simple equations based on the 
conservation of mass, energy, and momentum (cf. OC- 
CNP10): 

M acc = M inf - M outf , (6) 

where M acc is the mass rate actually accreted onto the 
BH. We define the kinetic energy rate of the outflow E w 

as, 

E w = e w M acc c 2 , (7a) 
= -M outl vl, (7b) 

p = M outf v w , (8) 

where we have oversimplified matters by allowi ng only 
one wind velocity, v w , when in fact Equation (|7b|) re- 
quires {v^) and Equation © requires (v w ) and e w de- 
notes the feedback efficiency, i.e., the wind efficiency in 
the momentum feedback model. Now, defining 

tp = 2e w c 2 /vl = M outf /M acc , (9) 
we have, as solutions to Equations ^ - flSJ), 

M acc = M inf — — , (fOa) 
1 + tp 

Moutf = Minf—--- , (10b) 
1 + tp 

E w = e w c 2 M inf —!—, (fOc) 
1 + tp 

P = Minf«w— 7T . (lOd) 
1 + tp 

As we see from Equations (fTOk) and (|10b). there is 
an important dimensionless quantity tp = M out f/M acc . 
In typical treatments of AGN feedback, tp is assumed 
to be in equations (fTU)) implicitly assuming v w —> oo, 

and so JW ou tf and p are neglected. Thus, the two terms 
E w and M acc may be overestimated, and p, the momen- 
tum input to the surrounding fluid is neglected. For ex- 
ample, if we adopt for the efficiency of generating me- 
chanical en ergy the value e w = e/e r = 5 x 10~ 3 , as 
adopted by iSpringell (|2005() and other authors, and we 
take i> w = 10, 00 km s" 1 (v w , w ) (jCrenshaw et all 120031 : 
iMoe et al.|[2009l ). the n we have tp — 9 v w 2 10 and all of 
the neglected effects may in fact be dominant; 90 %, 
{tp/(l + tp)) of the inflowing mass may be ejected in a 
disk broad absorption line (BAL) wind and the mass 
and momentum input deposited in the ambient gas may 
dominate over the energy input. 

To implement the actual output of mass, momentum, 
and energy, we modify the stochastic approach applied 
for the gas particle accretion on the BH shown in Equa- 
tion ([3]). We first calculate a probability of being at- 
tracted into the central zone by the BH for each gas par- 
ticle nearby using Equation Q and determine its fate 
by generating a random number Xj in the interval [0,1]. 
For Xj < pj, the gas particle is taken to be part of the 
inflow onto the BH. We then draw an independent ran- 
dom number tjj uniformly distributed in the interval [0,1] 



and compare it with 1/(1 + tp), the probability of being 
actually absorbed by the BH. For tjj < 1/(1 + tp), the 
gas particle is accreted onto the BH. For tjj > l/(l + tp), 
the gas particle is ejected with its wind velocity v w and 
momentum. 

We fix the wind velocity f w to 10,00 km s _1 (w w ,io) 
(|Crenshaw et al.ll2003trMoe et al.|[2009T) corresponding to 
a typical broad line wind. Together with our choice of 
e w = 5 x f 0~ 4 (note that our energy coupling parame- 
ter e w is 10 times smaller than the value used in other 
GADGET-2-based simulations), we have tp = 0.9. That 
is, essentially one of the two particles inflowing to the 
BH is actually accreted by it, the other is driven out as 
part of the broad absorption line wind. 

We emit a particle radially from the BH in the specified 
direction. We set the direction of wind to be parallel or 
anti-parallel to the direction of angular momentum (fx v) 
of each gas particle, making it to be essentially perpen- 
dicular t o the disk. The outflow i s stronger perpendicular 
to disk (jProga fc Kallm an 2004) as the feedback is rela- 
tively inefficient in the accreting disk of BH which sup- 
plies a continuous fuel to the BH via inflow. In addition, 
the massive and geometrically thick nature of the molec- 
ular torus (e.g.. iKrolik fc Begelmani ri988: Ta cconi et al.l 
I1994T ) may restrict the outflow to the direction essentially 
parallel to the disk angular momentum vector. 

AGN outflows can collide with and shock ambient gas, 
generating a momentum-driven flow. To mimic this phe- 
nomenon, we let the emitted wind particle share its mo- 
mentum with the neighboring gas particles. With this 
"momentum share" treatment, the two nearest neigh- 
boring gas particles are expelled together with the wind 
particle. They have the same velocity increment, Au w ~ 
10,000/3 km s _1 , conserving the momentum. Sharing 
momentum with other particles via inelastic collisions, 
however, d ecre ases the total kinetic energy increment 
(Equation J7b|) while preserving momentum. We de- 
posit the residual energy into these three particles in 
thermal form so that total energy is conserved. Note 
that the wind particles can reach very high temperatures. 
The analysis of the number of momentum sharing parti- 
cles will be discussed in the later sections. 

Momentum sharing has two technical advantages. In 
general, high velocity particles will drive shocks. Sim- 
ilarly to SN remnants, the solution will approach a 
Sedov solution. However, with coarse mass and time 
resolution and having momentum share starts the cas- 
cade of transforming the ~ 100% kinetic energy flux 
of the wind outflow to the ~ 25% kinetic energy flux 
of the Sedov solution with twice the number of parti- 
cles, a lower initial fraction of kinetic energy; and it 
makes the approach to the Sedov solution faster. Thus 
it makes us less sensitive to the problem of not hav- 
ing enough particles to correctly represent a hydrody- 
namic outflow. It also has computational advantages 
with regard to the time stepping. The standard Courant 
time step calculation implemented in the public release 
of GADGET might not guarantee fin e enough tim e 
stepping for strong explosion problems (|Springeli r2010'). 
and requires additional t ime-step limiter implem entation 
(jSaitoh fc Makinoll^OOHrDurier fc Vecchiall2012ft . In mo- 
mentum sharing, reducing the velocity of the wind by a 
factor of three, while maintaining the same momentum 
flux, reduces the need for short and expensive time steps 
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by the same factor, and our experiments have shown that 
if we perform solutions with high enough time and mass 
resolutions there is almost no difference between the two 
models. However, there is a great saving computation- 
ally in doing a momentum sharing when we are at low 
resolution. 

In this approach mass, energy, and momentum are 
transferred to the ambient gas during AGN feedback and 
in addition to the usual efficiency parameter e w , we must 
introduce a parameter ip (Equation ((9J), the ratio of the 
outflow mass flux to the accreted mass flux, which is fixed 
by e w and the wind velocity v w . 

2.2.3. X-ray Feedback 

We also consider the radiative feedback from the elec- 
tromagnetic energy component of X-ray radiation from 
the BH. We first calculate the total radiation emitted 
from the location of the BH particle as 



L r = e r M B HC , 



(11) 



where the ra diative efficiency e r = . 1 is adopted in all 
simulations dCiotti fc Ostrikerl 120071 ; iCiotti et al.l 120091: 
iNobel et al.ll2009l) . This luminosity is converted to a lu- 
minosity flux at the position of each particle by F r = 
L r /inr ■ , where r denotes the distance of the particle 
from the BH particle. We convert the flux to the net vol- 
ume heating rate E by using the iSazonov et al.l (2005) 
formulae, which describe the net heating rate per unit 
volume of a cosmic plasma in photoionization equilib- 
rium with a radiation field characterized by the aver- 
age quasar spectral energ y distribution, as in lCiotti et al.l 
(|2010fklNovak et al.l (|2011l ). We take into account Comp- 
ton heating and photoionization heating as summarized 
below. The volume heating rate E in cgs units is given 
as: 



E = n 2 {S 1 +S 2 ), 



(12) 



where n is the proton density (in number). The Compton 
heating term S\ is 



Si = 4.1 x 10~ 35 (1.9 x 10 7 - T) f , 
where the ionization parameter £ is defined as 



n(r)r 2 

The sum of photoionization heating is given as 

& = 10 - 23 a m ) b 



where 



6 = 1.1- 



i + (£/£o) b ' 

1.7 x 10 4 

J>0.7 ' 

1.1 4xl0 15 



e T/1.810 5 



T 4 



(13) 
(14) 

(15) 

(16) 
(17) 



and finally 



6> = 



1.5/VT+ 1.5 x 10 12 /\/T5 



4 x 10 1 

J>2 



Ml 



o (T-10 4 )/1.5 10 3 



(18) 



A speed-of-light delay in propagation from the BH is 
not included, since the effects of the delay should be small 
because of the small simulation scale (< 50 kpc). We 
neither include radiative transfer nor optical depth effects 
for the hard X-ray radiation considered here. 

However, we do include the electromagnetic 
momentum — the radiatio n pressure fr o m th e X-ra y 
flux from the BH as in iDeBuhr et ail 1)20101 l2011al) . 
We model the radiation pressure by applying a total 
momentum per unit time of 



P = E/c, 



(19) 



away from the BH particle to each gas particle, where 
E is the ene r gy ab sorbed by the particle given the 
ISazonov et al.l (|2005f ) prescriptions. Here we neglect the 
effect of dust since the IS M generally has a low opti- 
cal depth t o har d X-rays. lHambrick et al.l (|2011l ) and 
IKim et all (|2011D recently studied the effects of X-ray 
radiation on the properties of massive elliptical galaxies. 
They found that X-ray feedback is more effective at sup- 
pressing star formation and BH mass growth compared 
to the traditional thermal feedback model. 

2.3. Eddington Force 

In the previous studies, the maximum accretion has 
been limited to the Eddington rate (Equation ((2J) as 
shown in Section 12.21 Instead of manually limiting the 
mass accretion, we actually compute the Eddington force 
(EF) acting on the surrounding gas particles, directed 
radially away from the SMBH. We first calculate the lu- 
minosity as in Equation (jlip , and the flux at the position 
of the each particle by F r = L r /4irr 2 where r denotes the 
distance of the particle from the BH particle. Then, the 
total momentum change per unit time by the EF act- 
ing on the gas particles radially away from the SMBH 
particle is given as 



P : 



F r N e a T 



(20) 



where N e denotes the number of electrons associated 
with each gas particle and o~t is the Thomson cross- 
section for the electron. When the SMBH has vigorous 
mass accretion bursts, i.e., above the Eddington mass ac- 
cretion limit, the strong radiation pressure by the SMBH 
pushes the gas particles away, resulting in a lower den- 
sity near the SMBH particle, which leads to the lower 
mass accretion rate. Since the Thompson cross-section 
is independent of frequency (hv -C m e c 2 ), no radiative 
transfer treatment is required; we assume that absorbed 
UV radiation is re-radiated as IR radiation. 

2.4. New BH accretion model 

2.4.1. Bondi Radius Criterion 

The key to understanding the accretion process lies in 
correctly modeling the behavior of the accreting gas once 
it falls within the gravitational influence of the BH, the 
Bondi radius, reondi, defined as 



^Bondi 



2GM BH 



(21) 



where v 2 = c 2 + w rc i 2 , c and i> rc i denote the sound speed 
and the relative velocity to the SMBH of the gas respec- 
tively. The Bondi radius divides the flow into two distinct 
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regions. Far outside of rBondi, gas is hardly aware of the 
existence of the BH, and the flow is subsonic. On the 
other hand, inside the Bondi radius, gas has negative to- 
tal energy and essentially plunges at free-fall. The Bondi 
radius in the spherical case being the place where the 
Mach number of the flow is unity. We note that heating 
that occurs within the Bondi radius cannot, by defini- 
tion, affect the accretion rate since information cannot 
propagate upstream in a supersonic flow. 

Applying the Bondi-Hoyle-Lyttleton formalism 
(Equation (p}) to the BH mass accretion rate assumes 
that the accretion onto the SMBH is commensurate 
with the accretion rate through the Bondi radius. In 
cases where the SPH smoothing length is greater than 
the Bondi radius, i.e., the inner flow is numerically un- 
resolved, unusual effects can occur. To see what would 
happen if we used the standard accretion algorithm 
without feedback, we have made some artificial tests 
without any feedback mechanisms. We obtain a very 
large growth of the SMBH with M BH ~ 10 10 M Q . The 
neighboring gas particles around a BH at any distance 
from the BH are accreted, even for particles which are 
not gravitationally bound to the BH. 

To avoid any unphysical accretion from outside of the 
Bondi radius, we limit the accretion of mass to the gas 
particles statistically within the Bondi radius. With this 
Bondi radius criterion, gas particles can only be accreted 
when n < ^Bondi, where is the distance of the gas 
particle from the BH particle. When the mass of the 
BH is small, we cannot resolve the Bondi radius, i.e., 
the smallest resolvable length scale of our simulations, 
the gravitational softening length, is of a few tens of pc, 
whereas the Bondi radius is just a few pc when the BH 
mass is around 10 5 — 10 6 M Q (reondi/pc = 3.4 Mbh, &/v 2 50 
where Mbh,6 = 10 6 Af Q and v i50 =50 km s" 1 ). In this 
case, we set the Bondi radius to be the smallest resolved 
scale, i.e., the gravitational softening length of the gas 
particles. 

Since the gas mass distribution is smoothed with a ker- 
nel size hi, we apply what we term "a soft cut of Bondi 
radius criterion" . We allow the full accretion rate, only 
when the distance of the gas particle from the SMBH 
particles r, is smaller than rBondi + r so ft- The gravi- 
tational softening length r so ft, which avoids numerical 
singularities in the integral representation of the poten- 
tial, is the smallest resolvable length scale and serves as 
the minimum bound to the smoothing length hi. The 
soft Bondi probability psb as a function of the particle 
distance from the BH, for two cases of rs < ''soft and 
rs > r so ft, is shown in Figured] When the soft Bondi 
limit (SB) is included, we reduce the probability of be- 
ing absorbed by the BH for each gas particle (given as 
Equation ([3])) used in the original code by the soft Bondi 
probability psb, i-e., the final probability will be given 
as penal = Pj x Psb- This treatment essentially limits 
the mass accretion only to the gas particles statistically 
within the Bondi radius. 

2.4.2. Free-fall modification of accretion rate 

In the standard version of the code the actual accre- 
tion of the gas particles is determined by the probability, 
which is only a function of the SPH smoothing kernel. 
We have altered this to include a dependence on the time 
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Fig. 1. — Soft Bondi probability of being absorbed by BH as a 
function of the particle distance r. The probability for the case 
when the Bondi radius rBondi is equal to the softening length r so f t 
shown in the upper panel, and one for the case when rBondi is 
larger than r so f t is shown below. 

that it would take for the particle to be accreted allowing 
an extra factor of 



Pj,e 



(22) 



— T - 



where Tj = rj/(c 2 + v 2 ) 1 / 2 is the free fall time and N sp ^ 
denotes the typical number of smoothing neighboring gas 
particles of the BH. We modify the probabilities (Equa- 
tion to make them proportional to pjfi , i.e., we make 
it more likely that nearby particles will be accreted in 
a given time step than more distant ones. When the 
free-fall (FF) modification is included, the final proba- 
bilities of the particles to be accreted by BH are given as 
Pflnai = Pj x Pj, fi- where pj denotes a probability of being 
absorbed by BH used in the original code (Equation (J3J). 

2.4.3. Alternative averaging in the mass accretion 
calculation 



In the original implementation of iSpringel (120051) . the 
BH mass accretion is calculated based on the physical 
quantities such as density, sound speed, and relative ve- 
locity of the surrounding gas of the BH as shown in Equa- 
tion ([T]). These physical quantities for the BH particle 
are calculated from the ~ 64 neighboring gas particles 
using the SPH kernel. We can rewrite the Equation ([1]) 
as, 



Mr = 



47raG 2 M| H (p) 



)2)3/2 ' 



(23) 



where angle brackets denote the averaging over the neigh- 
boring particles using SPH kernel. This method can be 
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unstable to the number of SPH neighboring particles cho- 
sen and it separately averages quantities in the numera- 
tor and the denominator of Equation (I23|) . Suppose we 
made a run which used the 128 nearest particles instead 
of 64 nearest particles. It would have a different evo- 
lution, the soundspeed would be lower and so would the 
density. As noted earlier it is not the number of particles, 
per se that matters. If we halved the mass resolution and 
then moved to 32 particles we would keep the total mass 
into which the feedback energy was deposited a constant. 
We have performed this experiment and find that results 
are essentially unchanged. However, of course putting 
the energy into increasingly more mass would lower the 
feedback induced increase in Cg and change the Bondi 
rate after the use of Equation ([l) . 

Our new averaging method for the calculation of the 
BH mass accretion "alternative averaging (AA)" does 
the calculation in both time and space on an individ- 
ual particle basis and then averages the results over the 
nearest 64 particles in order to reduce the dependency on 
the number of SPH particles. With AA, we can rewrite 
Equation ([1]) as, 

m b,aa = ^ {ci+v2)m (24) 

where angle brackets denote the averaging over the SPH 
kernel. This method of performing the averaging is con- 
vergent because the added outer particles add progres- 
sively less and less. 

2.4.4. Fiducial BH mass accretion model 

In the fiducial model, we include all of the modifica- 
tions we described for the BH mass accretion model, i.e., 
soft Bondi radius criterion (SB), free-fall modification 
(FF), and alternative averaging (AA). We first calculate 
the BH mass accretion rate using Equation (|24p as shown 
above, and for the actual accretion of the gas particles 
onto the BH particle, we calculate a probability of being 
absorbed for each gas particle j as 

Pj, final = Pj X p jlS B X Pj,S- (25) 

We test our modified BH accretion model against the 
analytic Bondi solution in an idealized homogenous envi- 
ronment. As demonstrated in the Appendix we recover 
the Bondi solution. 

3. RESULTS 

3.1. Galaxy Models 

The disk galaxies used in our study are set up in dy- 
namical equilibrium and consist of a dark matter halo, a 
rotationally supported exponential disk of gas and stars, 
and a central bulge. The details of the model construc- 
tion are summarized in SdMH05. We test and study 
the stability of the constructed galaxy feedback models 
with a representative sample of our disk galaxy mod- 
els in isolation with v v - lr — 160 km s _1 , and r v ; r = 160 
hr 1 kpc corresponding to a virial mass of M v - lT = 9.53 x 
10 11 Ii~ 1 Mq. The dimensionless Hubble parameter is 
h = 0.71 such that the present- day Hubbl e param eter is 
Ho = 71 km s" 1 Mpc- 1 . The iHernquistl (fl990l) profile 
dark matter halos are constructed with the concentration 
parameter c = 9 of the corresponding Navarro-Frenk- 



White (NRW) halo (|Navarro et al.lll997l ). The dark mat- 
ter halo is then populated with exponential disks with a 
baryonic mass fraction of m,d = 0.041, so that a total 
disk mass M<j = mdM v - lv with a fractional gas content of 
/gas = 0.2 with the rest being stars. 

To study the effects of numerical resolution, we use 
four models with different mass and spatial resolutions 
but with the same initial setup. The resolution details, 
including the number of particles, particle mass, and 
corresponding gravitational softening lengths are given 
in Table 1. The dimensionless accretion parameter in 
Equation (P) was set to a = 100 in the previous stud- 
ies (SdMH05; Uohansson et a l. 2009b), which correspond 
to our low-resolution model. This value is quite a bit 
higher than the theoretical value of a ~ 1. This dis- 
crepancy is due to the numerical resolution limits, that 
is, the gas density and sound speed at the location of 
the SMBH are estimated over the large-scale, resulting 
in artificially low densities and high sound speed. The 
higher value of a has been adopted empirically to correct 
for this resolution effect. For the model with the higher 
numerical resolution, we adopt the smaller value of the 
dimensionless accretion parameter that can result in the 
similar scale of early accretion history of low resolution 
model with a = 100. We run a series of simulations with- 
out any feedback prescription with the different values of 
a and adopt the a value that best reproduces the early 
accretion history of the low-resolution run with a = 100. 
We turn off all the BH feedback to remove the resolution 
dependency of the feedback prescription. The adopted 
values of a for each resolution are shown in Table 1. 

We summarize the properties of the models in Table 2. 
For reference, we first list the thermal feedback model 
"T" , with the standard BH mass accretion and feedback 
model described in SdMH05. Note that model properties 
of standard model "T" in this paper are different from 
those of SdMH05, as we adopt "Very High" resolution as 
a standard resolution in this study, and the energy cou- 
pling efficiency e w = 5 x 10~ 4 which is 10 times smaller 
than the value adopted in SdMH05. We name our best 
proposed model with all the modifications that we have 
described "Fiducial" . To isolate and compare the effects 
of each feedback model or modification, we include some 
simulations with one modification missing. We also in- 
clude the "No Feedback" model: with no BH feedback in 
any form. Details of each model are shown in Table 2. 

3.2. Comparison of the different feedback mechanisms 

We first examine the effects of various feedback mecha- 
nisms on the evolution of the BHs and galaxies. Figure^] 
shows the global star formation rate, the accretion rate 
onto the SMBH, and the evolution of the BH mass for 
the two different feedback models, the classical, thermal 
feedback (e.g., SdMH05), and our best proposed "Fidu- 
cial" model which has mechanical feedback with the X- 
ray heating and radiation pressure. Note that both mod- 
els adopt the same feedback energy coupling efficiency 
e w = 5 x 10~ 4 , which is 10 times smaller than the value 
adopted in SdMH05. A seed SMBH starts at rest in 
the center with mass of Mbh = 10 5 Af© in all models, 
and grows due to gas accretion during the simulation. 
Note that the growth of BH is essentially the same in 
the two models. Thus the primary results of SdMHOS 
and lDi Matteo et "ail (|2005|) with respect to the growth of 
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TABLE l 

Summary of the Numerical Resolution 



Model 


DM Particles 


Disk Particles 


m DM {h~ 1 M Q ) 




a 

fc bar 


e DM 


a b 


Low c 


30,000 


40,000 


2.96 x 10 Y 


3.91 x 10 5 


0.1 


0.8 


100 


High 


400,000 


300,000 


2.25 x 10 6 


1.30 x 10 5 


0.02 


0.083 


35 


Very high 


800,000 


600,000 


1.13 x 10 6 


6.50 x 10 4 


0.016 


0.066 


32.5 


Ultra high 


1,600,000 


1,200,000 


5.62 x 10 s 


3.25 x 10 4 


0.013 


0.052 


30 



a Gravitational softening lengths in h 1 kpc. 

b Dimensionless accretion par ameter in Equat ion JTJ. 

c Numerical resolution used in Springcl ct al. (2005b) corresponds to our low-resolution model. 



TABLE 2 
Summary of Model Properties 







Feedback 












log 


log 


log 


log 






Model 


T/M 


X-Ray 


X-Ray RP 


FF 


AA 


SB 


EF 


NEL 


AA/bh 


AM wind 


/eff a 

'bh 


T ■ _, b 

-'-'wind 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


0) 


(10) 


(11) 


(12) 


(13) 


(14) 


(15) 


1 


Thermal d 


Thermal 
















7.47 


8.92 


-2.80 


36.99 


53.17 


2 


Fiducial* 3 


Momentum 


V 


V 


V 


V 


V 


s/ 


V 


7.38 


8.24 


-2.66 


40.72 


1828.1 


3 


Fiducial w/o XRP 


Momentum 






V 


V 


V 


V 


V 


7.96 


9.22 


-2.22 


41.37 


1287.9 


4 


Fiducial w/o EF 


Momentum 


V 


V 


V 


V 


V 




V 


7.41 


8.56 


-2.51 


41.07 


1727.2 


5 


Fiducial w/o SB 


Momentum 


V 


V 


V 


V 




V 


V 


7.54 


8.71 


-2.37 


41.08 


1743.9 


6 


N 


No feedback 
















9.81 




-3.86 






7 


N-SB 


No feedback 










V 






8.89 




-2.38 







Notes. Model names indicate the activated physics (symbol ^/) in the simulations as detailed in Column 3—10. For example, in Model 
"Thermal" only thermal energy feedback is allowed, while in Model "N-SB", no feedback is included, with the soft Bondi (SB) treatment. 

a 'l3H — ^bh opt/^Edd where i B 1j opt is the BH luminosity in the optical band after absorption, i.e., as it will be seen from infinity. The 
mean Eddington rates for last 1 Gyr are listed. 

b The mechanical luminosity of the wind based on the wind mass rate and wind velocity measured at 5 kpc from BH for last ~ 1 Gyr. 

iwind = Af wind ^ ind /2. 

c Gas wind velocity at 5 kpc from BH in km s . 

d This model corresponds to the purely thermal feedback model discussed in SdMH05 (Springcl ct al. 2005b). 
e Our best proposed "Fiducial" model. 



central BHs during galaxy mergers are not substantially 
altered by the changes that we have introduced. Also, 
overall star formation rates of the two models are simi- 
lar. However, the new model with mechanical feedback 
has more episodic star formation as in the mass accre- 
tion rate. Accretion rates and radiation output are much 
more variable in the new treatment with episodes of high 
accretion now reaching ^BH/^Edd > 0.1. The mechani- 
cal feedback model spends a large fraction of time at rel- 
atively low Eddington accretion rate, coinciding with the 
observed broad-line active galaxies in the local universe 
(iGreene fc Hoi 120071 : iHol 120091 : iKauffmann fc Heckmanl 
2009). 

Moreover, the two models have quite different wind 
properties. In the mechanical feedback model, wind par- 
ticles are ejected with the instantaneous disk wind veloc- 
ity of v w ~ 3000 km s" 1 (OCCNP10), while the heating 
from AGN feedback energy drives slow and hot outflows 
from galaxies in the thermal feedback model. The exis- 
tence of a weak wind perpendicular to the plane of the 
disk in the vicinity of the BH was shown in SdMH05. 

In order to compare the wind properties, we first mea- 
sure the total outflow wind mass throughout the sim- 



ulation by summing up the mass of the gas particles 
which pass through the plane \z\ = 5 kpc, above and 
below the disk midplane. Then, we measure the ve- 
locity of the particle at each time step and calculate 
the corresponding mechanical luminosity as £ w i n d = 
-^wind^ ind /2, i.e., the kinetic energy carried away by 
the outflowing winds. Temporal evolution of the wind 
properties (i.e., outflow rate, wind velocity and mechan- 
ical luminosity) of two representative models of ther- 
mal feedback (SdMH05) and the fiducial model is shown 
in Figure [3] Once the BHs reach similar masses, af- 
ter about 3 Gyrs, the thermal feedback model devel- 
ops a wind with about a factor of 10 higher outflow 



rate, but the velocity of the wind v v 



100 km s" 



is a factor of 20 smaller compared to the mechani- 
cal feedback model. The high wind velocity (v w ~ 
2000 km s _1 ) in our fiducial model is consistent with the 
recent v elocity measurement of mass outflows in local 
Seyferts (jFischer et alJl20lH I Miiller- Sanchez et~alll20lH 
iPounds fc Vaughan 2011) which ranges from 700 km s _1 
up to 3000 kms" 1 . 

The mechanical luminosities of the winds for two mod- 
els are shown in the middle panel of Figure [U and 
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Fig. 2. — Comparison of the feedback models: classical thermal 
feedback (blue), and our best proposed "Fiducial" model which 
has momentum mechanical feedback with the X-ray heating and 
radiation pressure (red). Evolution of the global star formation 
rate (top), the accretion rate onto the black hole (middle) and the 
evolution of the black hole mass (bottom) for an isolated gas-rich 
disk galaxy. Note that both the star formation rate and the black 
hole growth are essentially similar in the new momentum driven 
treatment and the previous thermal feedback model. However, the 
degree of fluctuation in the mass accretion rate is greater in the 
momentum feedback model. 

the averaged mechanical luminosities for last 1 Gyr 
are listed in Table 2. The thermal feedback model 
has considerably smaller mechanical luminosity L w - m ^ ~ 
10 37 erg s _1 , because of its slow wind velocity. On the 
other hand, our fiducial model with mechanical feedback 
has strong outflow with L W i n d ~ 10 41 erg s -1 , which 
is consistent with the fact that large outflows with a 
kinetic power corresponding to a significant fraction of 
the AGN bolometric luminosity are commonly observed 
in X-ray observations of a number of quasars (mostly 
BAL systems) that reveal significant absorption (e.g., 
Chartas et aTIl2003HCrenshaw et all 120031: iPounds et al.1 
2003tlHolt et al.ll2008[ ). The total kinetic energy carried 
away by the winds, i.e., the mechanical luminosities in- 
tegrated over the entire simulation time for the thermal 
model is Ai? w ; n d ~ 2.8 x 10 55 erg, while the momentum 
feedback model deposits Ai? win d ~ 8.3 x 10 57 erg into 
the ISM within 6 Gyr. 

From a feedback energy coupling efficiency e w which 
is an input parameter, and a total mass change for the 
BH, we calculate the energy released by the BH for the 
two cases using A£J mec h = CwAA/bhc 2 and compare the 
ratios of the energy released in winds to the energy re- 
leased by the BH, A£' w i n d/A£ , mcc h. While the thermal 
feedback model emits only 0.1 % of the BH mechanical 
energy release in winds (AE win( i / AE mcc h = 0.0011), the 
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Fig. 3. — Top: the evolution of wind mass loss rate; middle: cor- 
responding mechanical luminosities; bottom: gas wind velocities. 
All quantities are measured at 5 kpc from the SMBH, above and 
below the disk midplane. Note that the thermal feedback model 
has much smaller wind velocity, a factor of 20 smaller compared to 
the mechanical feedback model. 

value rises to 38 % for the momentum feedback model, 
Ai? w i n d/Ai? m ech = 0.38. That is, the mechanical energy 
put into momentum drives winds more efficiently than 
the energy put into heat. 

Now we investigate the behavior of the two models 
for disk galaxy simulations at higher numerical resolu- 
tion. However, such a study is complicated because of 
numerical resolution issues. In the previous studies using 
one-dimensional simulations, it was found that including 
the mass and momentum component has the largest ef- 
fect on the final mass of the SMBH, reducing the final 
SMBH mass by a factor of up to 100 (OCCNP10). Turn- 
ing on or off the energy input has relatively small effect, 
altering the SMBH mass only by a factor of two. But, 
in this three-dimensional study we have found that the 
momentum feedback is not more efficient than thermal 
feedback in protecting the SMBH from growth. However, 
the three-dimensional classic treatment has resolution ef- 
fects that are difficult to correct because (1) the feedback 
energy is deposited outside of the Bondi radius whereas 
it is distributed within the Bondi radius in the much 
higher resolution one- and two- dimensional studies, and 
(2) the accretion is determined with the estimated gas 
density and sound speed averaged over the smoothing 
kernel, which is perhaps not the optimal procedure. 

The results of the resolution dependency and the com- 
parison of the wind properties of two feedback models are 
shown in Figure @] Note that we adopt the smaller value 
of dimensionless accretion parameter a for the higher 
resolution runs as listed in Table 1. As described above, 
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Fig. 4. — Outflow wind properties and the final SMBH mass for 
the thermal feedback (blue) and the momentum feedback (red) are 
shown as a function of resolution (gas particle mass). The total 
outflow wind mass throughout the simulations, the averaged out- 
ward wind velocity, the averaged mechanical luminosity for last 1 
Gyr, and the final SMBH mass (from top to bottom) are shown. In 
the thermal feedback model, much weaken winds are generated in 
the highest resolution, since the injected feedback energy is quickly 
radiated away due to the very short cooling time of the gas. 

we measure the outflow wind properties throughout the 
simulations at \z\ — 5 kpc from BH, above and below 
the disk midplane. The measured wind masses, the time 
averaged wind velocities, and the averaged mechanical 
luminosities over the last 1 Gyr are shown for each feed- 
back model as a function of the gas particle mass in Fig- 
ure UJ The total amount of mass in the outflowing wind 
is resolution dependent in both models, but the effects 
of resolution appear to be larger in the thermal feed- 
back model. While the momentum feedback model in 
the highest resolution simulation has less wind mass (by 
a factor of 10) compared to the lowest resolution one, the 
wind mass difference between thermal feedback models 
at different resolutions reaches 10 3 . The wind velocity is 
resolution-independent for both feedback models, how- 
ever thermal feedback models have much slower winds 
(~ 50 — 100 km s _1 ) compared to momentum feedback 
in all resolutions. The AGN-driven wind in the ther- 
mal feedback is also much slower than the velocity of the 
observed broad absorption line winds ~ 10,000 km s _1 
(jCrenshaw et al J 1200a iMoe et all [20091 ) . In the case of 
the thermal feedback model at higher resolution, the 
shorter cooling time of dense gas makes thermal input 



increasingly inefficient at higher and higher numerical 
resolution. 

In the case of the final mass of the BH (bottom panel of 
Figure 0J, the effect of the resolution in the momentum 
feedback model is moderate, whereas the thermal feed- 
back model has a factor of 10 smaller final mass in the 
highest resolution compared to one in the lowest resolu- 
tion. For pure thermal feedback, we deposit the thermal 
feedback energy into the neighboring gas particles of the 
BH. The number of the affected neighboring gas parti- 
cles is set to be the same for the different resolutions, 
therefore, in the higher resolution studies, we add energy 
into a smaller mass of the gas. In the higher resolu- 
tion cases, gas particles in the central region have higher 
sound speed and lower density resulting in the lower BH 
accretion rate. However, the thermal feedback depends 
on the mass into which the energy is deposited, not on 
the number of particles. If the feedback energy were 
spread over constant mass by increasing the number of 
particles into which the thermal energy were added in the 
higher resolution study, the results would be essentially 
the same. 

We compare the evolution of hot gas a nd X-ray emis- 
sion o f the two feedback models. Following lNavarro et al.l 
(1995ft . we assume that X-rays are produced by the cool- 
ing of hot and diffuse gas, and estimate the X-ray lumi- 
nosity for each SPH gas particle as, 



I Xi = 1.2 x 10 



_24.Pi^gas,i ( ^-^i \ 

(nm p ) 2 V IkeV / 



1/2 



erg s 



(26) 



where pi, TO gaSj i and Ti are the density, mass, and tem- 
perature of the i th gas particle in cgs units, respectively, 
m p is the proton mass, and p is the mean molecular 
weight. The X-ray emission computed via Equation (|26|) 
for each SPH particle is a lower limit as it assumes that 
the primary mechanism for X-ray emission is thermal 
bremsstrahlung. We do not include the X-ray emission 
via metal-line cooling although it is the more efficient 
cooling mechanism for metal-enriched gas with a temper- 
ature of ~ 10 6 K and would produce more X-ray emis- 
sion. We also assume that the central region of the galaxy 
remains obscured because of the large column density of 
intervening gas and dust, and calculate the total X-ray 
luminosity as 



iv h 



-'X.tot 



E 



J X,i, 



(27) 



where th e sum is over all hot and diffuse gas particles. 
Following lCox et al.l (|2006h , we define the 'hot and diffuse 
gas' with a temperature of T > 10 6 K, and a density p < 
3.16 x 10 _3 M Q pc~ 3 , which corresponds to the critical 
density for star formation. 

In Figure O we show where the SPH particles lie in 
the density-temperature plane with the calculated X-ray 
luminosities for both feedback models. In the thermal 
feedback model, the energy input by the accreting BHs 
generates collimated, hot and slow winds perpendicular 
to the disk plane. These wind particles leave from the 
high density and hot temperature tip of the distribu- 
tion of gas in the density-temperature plane, and cool 
down slowly as they move outward. Due to their low 
velocity and high density, the cooling time is relatively 
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Fig. 5. — Distribution of gas in the density-temperature plane at t = 5.95 Gyr for the two feedback models. The radial distance from 
the BH is color coded, and the number density of SPH particles on the plot is shown with contours. The horizontal and vertical dotted 
lines respectively represent the temperature and density cuts we adopted for the X-ray luminosity estimate. Note that compared to the 
momentum feedback model, the thermal feedback model has much higher X-ray luminosity Lx ~ 3 X 10 4t) erg s — 1 , which is emitted from 
the diffuse hot gas, even during the low BH accretion phase (LBH/^Edd ~ 0.002). 



long, and these hot and diffuse winds emit X-rays with 
the luminosity far greater ( ~ 10 40 erg s _1 ) than the 
momentum driven winds (~ 10 36 erg s _1 ). The mo- 
mentum feedback model has much lower X-ray luminos- 
ity as the momentum-driven winds quickly expel the X- 
ray emitting hot gas. Even during the low BH accre- 
tion phase with LBH/^Edd ~ 0.002, the X-ray luminos- 
ity of the thermal feedback model is much higher than 
what is typically seen from normal massive spiral galax- 
ies at the pre sent epoch (10 38 -10 39 erg s" 1 in the 0.5 



2.0 keV band (lAlmv et al. l l2000HOwen fc Warwick |[200l 



iBoroson et alJl2011l )). The X-ray luminosity of the ther- 
mal feedback model is a lower bound, since the inclusion 
of cosmologically infalling gas in the simulations would 
lead to an increase in the computed thermal X-ray emis- 
sion. 

We now discuss the effect of the amount of the 
momentum-driven flow, i.e., number of momentum shar- 
ing neighbors in the momentum feedback model. As de- 
scribed in Section 12.2.21 we let the wind particle share 
its momentum with the neighboring gas particles, to 
mimic the shocked swept-up ambient medium. Because 
of the resolution limit, we keep the number of momen- 
tum sharing neighbors constant, despite the fact that a 
total swept-mass by momentum- driven wind depends on 
the mass of the BH (jKingl 12003). In order to study the 
effect of the number of neighbors to share the momen- 
tum feedback, we ran a series of simulations adopting 
different numbers of neighbors. We find that the effects 
of sharing momentum with more gas particles are small. 
Adding more particles is equivalent to assuming an early 
and brief transition to the Sedov phase. 



3.3. Galaxy Merger Simulation 

In addition to isolated galaxies, we also performed an 
equal-mass galaxy merger simulation using our progen- 
itor disk galaxy models. The merger simulation was 
ran at high numerical resolution (see Table 1) with 
the initial seed BH masses set at 1O 5 M . Follow- 
ing Uohansson et al.l (l2009af) we adopt orbital geome- 
try G13 (|Naab fc Burkertl 12001 ) for our merger simu- 
lation. This geometry corresponds to the inclinations 
i p = — 109, « s = 180 and the arguments of the pericenter 
ujp = 60, u> s = for the primary and secondary galax- 
ies, respectively. The galaxies approach each other on a 
parabolic orbit where the initial separation of the pro- 
genitors is i?i n it = r v i r and the pericentric distance is 
r por i = 2r,i, where r v - lr = 160ft. -1 kpc is the virial radius 
and rd = 2.5/i -1 kpc is the disk scale radius. The sim- 
ulation was evolved for a total of t = 3 Gyr with the 
merger taking place at t ~ 1.5 Gyr. The equal-mass 
merger initial conditions are simulated using both the 
standard thermal feedback and the new momentum me- 
chanical feedback with the X-ray heating and radiation 
pressure. 

In Figure El we show the evolution of the resulting to- 
tal star formation (top), SMBH accretion (middle), and 
SMBH mass (bottom) for the two feedback models as 
a function of time. A similar evolution is seen in the 
star formation rate for both models, however the new 
feedback fiducial model has episodic outbursts of mass 
accretion as shown in the isolated galaxy model. Even 
though the two feedback models have different SMBH ac- 
cretion history, the final mass of black hole is essentially 
the same. Other aspects of the two simulations, such as 
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Fig. 6. — Comparison of the feedback models with a major 
merger of two galaxies: classical thermal feedback (blue), and our 
best proposed "Fiducial" model which has momentum mechanical 
feedback with the X-ray heating and radiation pressure (red). Evo- 
lution of the total star formation rate (top), the total accretion rate 
onto the black hole (middle) and the evolution of the black hole 
mass (bottom) are shown as a function of time for the 1:1 merger. 
The filled circles indicate the time of SMBH merger. Note that 
both the star formation rate and the black hole growth are essen- 
tially similar in the previous thermal feedback model and the new 
momentum driven treatment. However, the degree of fluctuation 
in the mass accretion rate is greater in the momentum feedback 
model as shown in the isolated galaxy case. 

the X-ray thermal luminosity, are significantly different, 
and this will be discussed in later papers. 

3.4. Soft Bondi criterion 

Next, we turn our attention to the Bondi radius cri- 
terion. In order to test the new soft Bondi mechanisms 
for limiting accretion when numerical resolution is less 
than optimal, we ran two artificial test runs without any 
feedback, i.e., "No Feedback" models, with and without 
soft Bondi criterion. As discussed earlier in the Sec- 
tion I2.4.H since any closest neighboring gas particles 
from the SMBH at any distance are considered as po- 
tentially accreting particles, gas particles keep accreting 
onto the SMBH and finally BH ends up swallowing all 
the gas in the host galaxy in the "No Feedback" model 
without soft Bondi criterion. 

With the Bondi radius criterion added, we limit the ac- 
cretion only to the gas particles statistically within the 
Bondi radius and we can prohibit the gas particles which 
are not within the gravitational influence of the BH from 
accreting onto the SMBH. In more realistic simulations 
which included feedback, the effects of our Bondi limita- 
tion are far less significant. 

In Figure [3 we show the resulting SFRs, BH accretion 
rates, and BH mass growth for these two models with and 



3.0 r 



7 2.5 r 




Time [Gyr] 

Fig. 7. — Comparison of the models with and without the Bondi 
radius criterion. Evolution of the global star formation rate (top), 
the accretion rate onto the black hole (middle) and the evolution of 
the black hole mass (bottom) are shown for each model. Note the 
order of magnitude reduction in the total accreted mass occasioned 
by including the Bondi limiter. 

without the Bondi radius criterion. Utilizing the Bondi 
radius criterion effectively regulates the mass accretion, 
reducing the final mass by a factor of 10. 

3.5. Fiducial model and the effects of other physics 

In order to isolate and compare the effects of each feed- 
back mechanism or modification made in the accretion, 
we perform the test simulations with one modification 
missing. In Figure |H1 we show the temporal evolution 
of SFR, mass accretion rate onto BH and BH mass of 
three distinctive models, the one without radiative heat- 
ing and radiation pressure (fiducial w/o XRP), the one 
without the Eddington force (fiducial w/o EF), and the 
one without the soft Bondi radius criterion (fiducial w/o 
SB), along with the fiducial model with all modifications 
as a reference. The radiative heating and radiation pres- 
sure is most effective among them, limiting the final mass 
of BH a factor of four. Compared to radiative feedback, 
both EF and Bondi criterion have minor effects on the 
final BH mass. 

4. DISCUSSION 

Absent a specific mechanism for transferring energy to 
particles surrounding an accreting BH, the AGN feed- 
back via mass, and momentum output have not been 
included in classic thermal feedback works in three- 
dimensional hydrodynamic simulations. In this study, 
we have included the momentum mechanical feedback in 
the SPH simulation code, GADGET-2. We also include 
a treatment of the feedback by the X-ray radiation which 
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Fig. 8. — Various model comparison. Evolution of the global 
star formation rate (top), the accretion rate onto the black hole 
(middle) and the evolution of the black hole mass (bottom) in 
simulation of an isolated gas-rich galaxy. The largest of the newly 
included effects is clearly due to the X-ray heating and radiation 
pressure (XRP). 



emanates from the BH and heats the surrounding gas in 
the host galaxies as well as radial momentum added to 
the fluid. We statistically limit accretion to gas parti- 
cles which are gravitationally bound to the central BH 
(Bondi radius criterion). 

A series of test simulations of isolated systems with 
new implementations of the momentum feedback and ra- 
diative feedback with new criteria on mass accretion are 
performed. Relevant quantitative properties of the mod- 
els are presented in Table 2, while the general results can 
be summarized as follows. 

1. Overall, the BH growth is quite similar in the two 
approaches, so the successful predic tion of the Mbh — o" 
relation by iDi Matteo et all |2005l . |2008) would be ex- 
pected to be maintained in the new approach. 

2. Our best proposed fiducial model with mechan- 
ical and radiative feedback by hard X-rays has much 
higher velocity outflows compared to the thermal feed- 
back model, with i> w ~ 2000 km s _1 and £ w ind ~ 
10 42 erg s _1 . The total emitted kinetic energy of mechan- 
ical feedback model is ~ 100 times higher than that of the 
thermal feedback model, even when the same feedback 
energy coupling efficiency is assumed and the BH growth 
is similar. The outflows found in our fiducial model have 



properties broadl y similar to those observed in some local 
Seyfert galaxies dAlatalo et al.|[20ill iFischer et al.M2011t 
IMuller-Sanchez et alJl20lHr ~ 

3. The hot gas produced by slow, dense, ther- 
mally driven winds emits an X-ray luminosity signifi- 
cantly greater (~ 10 40 erg s~ x ) than the momentum- 
driven winds (~ 10 36 erg s _1 ). This X-ray luminos- 
ity of the thermal feedback model is also far greater 
than what is typic ally seen from normal spiral galax- 
ies (~ 10 38 erg s- MAlmv et~aI]l2000HOwen fc Warwick! 
l2009h . 

4. In the mechanical feedback model, the fluctuation 
level in both radiant and wind outputs is considerably 
greater than in the standard thermal feedback model. 
While the thermal feedback model has a steady mass 
accretion with the Eddington ratio isH/^Edd = 0.001- 
0.01 throughout the simulations, the momentum feed- 
back model has stochastic bursts in the mass accretion 
with the Eddington ratio, which spans from ieH/^Edd = 
10~ 6 to 1. 

5. In an artificial model computed without any feed- 
back mechanisms (no feedback model), the SMBH grows 
to ~ 10 10 Mo accreting all the gas particles in the host 
galaxy. As noted, the standard prescription for accretion 
does not require the accreted particles be gravitationally 
bound to the central BH. However, the statistical imple- 
mentation of the Bondi radius criterion can effectively 
limit the accretion of the gas particles to gravitationally 
bound particles, reducing the final mass of BH by a fac- 
tor of 10. In more realistic models with feedback the 
differences are small. 

6. Radiative heating and radiation pressure on the 
ISM by photons emitted by the central BH moderately 
reduces the final mass of BH, by a factor of four. 

7. The growth of the BH is confirmed to be essentially 
the same in the thermal and momentum feedback models 
in an equal-mass galaxy merger simulations. 
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APPENDIX 
ACCRETION MODEL 

We test our modified model of unresolved accretion onto the BH summarized in Section 12.41 Our simulation set-up 
consists of a random distribution of gas particles in a periodic box with an accreting BH particle in the center. The 
size of the box is 1 kpc, and 20000 gas particles with a mass of 2750 M Q are used for the simulations. The mass 
of each gas particles is ~ 30 times smaller than that of our galaxy simulation model at standard resolution, and the 
gravitational softening length is ~ 5 pc. We choose Mbh = 4 x 10 6 Af Q and the corresponding initial Bondi radius 
(Equation (|21[) ) is reondi = 10 pc, with T gas =20000 K and Gaussian velocity distribution with v gas = 50 km/s and 
°Vgas = 5 km/s. No radiative cooling is considered and dimensionless accretion parameter a is set to be 1. 
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Fig. 9. — Evolution of the black hole masses (top) and m ass a ccretion rates onto BH (bottom) of the simulations of mass accretion models. 
The standard mass accretion model described in Section 12.21 in this paper (refer to SdMH05 for details), and new accretion model with 
all of modifications we described for the BH mass accretion model, i.e., soft Bondi radius criterion (SB), free-fall modification (FF), and 
alternative averaging (A A) (see Section |2.4> . The integrated new accretion model (red) agrees well with the analytic expectation despite 
the fact that the instantaneous accretion rate is often higher due to the intervals when the instantaneous rate is zero. 

We perform a direct comparison of our new accretion model with the SdMH05 model. In Figure [9j we plot the 
resulting BH accretion rate, and total BH mass for the simulation performed using the standard mass accretion model 
described in SdMH05, and compare it to the corresponding output of the accretion model presented in this paper 
(with FF, SB radius criterion, and AA following Equations (|24p and ([25])). We also plot the analytic solution of the 
mass a ccretion of a given physical properties, i.e., Bondi solution which is described in the analytical formula of Bondi 
(|1952D under the assumption of spherical symmetry and negligible angular momentum as, 

M B = A47rr| ondi/ 9 00 c Si00 , (Al) 

where t he dimension l ess pa rameter A depends only on the adiabatic index of the accreting ga s (for details, see [Bondi 
(|1952t) ; Uaniuk et all (|2009t )V For an assumed adiabatic index 7 = 5/3, A = 0.25 (|Bondilll952f ). 

The evolution of the BH mass as a function of time for the new accretion model agrees well with that of the analytic 
Bondi solution. The total accreted gas mass for the previously used accretion model is about two times larger. As 
shown in the bottom panels of Figure [9j in the standard mode the mass is accreted continuously whereas in our 
fiducial model we have discrete accretion events. That is mainly because the soft Bondi radius criterion prohibits the 
gas particles which are not within the gravitational influence of the BH from accreting onto the BH and stops the 
accretion before the gas particles are found statistically within the Bondi radius. The new model does not have mass 
accretion for the first ~ 0.2 Gyr, before the accretion flow is formed around the BH in the center. Even if the Bondi 
radius is not fully resolved as in many of the applications, we obtain the Bondi solution with our new BH statistical mass 
accretion prescription. Given the discrete nature of our particles the accretion rate is inevitably stochastic. However, 
the computed accretion rate closely follows the exact analytic solution when we adopt our statistical treatment. 



